Theoretical evidence for a dense fluid precursor to crystallization 
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We present classical density functional theory calculations of the free energy landscape for fluids 
below their triple point as a function of density and crystallinity. We find that for both a model 
globular protein and for a simple atomic fluid modeled with a Lennard- Jones interaction, it is 
free-energetically easier to crystallize by passing through a metastable dense fluid in accord with 
the Ostwald rule of stages but in contrast to the alternative of ordering and densifying at once as 
assumed in the classical picture of crystallization. 
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Crystallization is an intricate process of fundamental 
importance in many areas of physics, chemistry and en- 
gineering. The classical picture of crystallization from 
supersaturated solutions goes back to Gibbs and consists 
of the spontaneous formation of crystalline clusters which 
then either grow or shrink depending on the relative im- 
portance of the free energy gain due to the lower bulk free 
energy of the crystal cluster and the free energy penalty 
due to the surface tension between the two phases. In 
this picture, the local density is the only order param- 
eter: the crystalline cluster is (in general) denser than 
the fluid. In recent years, this picture has been called 
into question by simulation, theory and experiment for 
the particular and important case of the crystallization 
of globular proteins, ten Wolde and Frenkel (hereafter 
tWF) showed by means of simulation that the free en- 
ergy landscape of protein crystal clusters as a function of 
the number of atoms in the cluster and the " crystallinity 
" favored paths leading from no clusters to clusters with 
low order to ordered clusters over paths moving from no 
clusters directly to ordered clusters This picture was 
confirmed by Talanquer and Oxtoby^ and Shiryayev 
and Guntonj^] who showed using a parameterized van 
der Waals-type model of globular proteins that surface 
wetting did indeed lower the free energy of crystal clus- 
ters. More recently, the simple picture has also been chal- 
lenged by novel experimental investigations. Vekilov and 
co-workers have shown that, prior to crystallization, pro- 
tein solutions harbor metastable droplets of dense fluid 
and they have suggested that these droplets are necessary 
precursors of crystallization^, 0, |^ Q- The picture that 
emerges is one of the formation of metastable droplets of 
dense fluid which then subsequently crystallize. In this 
paper, we show by means of classical density functional 
theory calculations that there is an intrinsic free-energy 
advantage in first densifying into a metastable dense- 
fluid state and then crystallizing rather than following 
the classical path which goes directly from gas to crystal. 
Furthermore, our calculations suggest that a similar ad- 
vantage exists for fluids of small molecules, modeled here 
via the Lennard- Jones (LJ) interaction, thus indicating 



that this mechanism may underlie most crystallization 
processes. 

The starting point for our analysis is classical density 
functional theory (DFT) which is based on a theorem, 
due to Mermin, that the Helmholtz free energy of a clas- 
sical system is a unique functional F [p] of the local den- 
sity p(~f*). The local density is a constant, p{~r') = p , 
for a bulk liquid while for a simple bulk solid it is a sum 
of localized functions centered on the lattice sites, 

1=0 

for some function / (T^)where the vectors ^^'C the 

lattice vectors and the lattice constant is oq . Typically in 
a bulk solid, this is approximated as a Gaussian, / (T*) = 

?7o (^)"^^^exp (— ar^), where ryg G [0, 1] is the fraction of 
lattice sites which are occupied and the parameter a is 
related to the degree of crystallinity. The average density 
for a lattice with Nq lattice sites per unit cell is p = 
V Iv P ~ VoNoa^^ where V is the volume of the 

system. Given the Gaussian approximation the density 
can also be written in terms of Fourier components as 

P(^) =P + P^(^^P [f^i ■ 'r'/a-oj cxp{-Kf/4:ala) , 
1=1 

(2) 

where ^-re the reciprocal lattice vectors . This form 

shows clearly that as a goes to zero, the density becomes 
uniform corresponding to a fluid whereas the real-space 
form shows that as a goes to infinity, the density becomes 
infinitely localized as a sum of Dirac delta functions. For 
this reason, it is natural to take x = exp(— i^^/4a) to 
be an order parameter corresponding measuring "crys- 
tallinity " since it becomes zero for the liquid and one for 
the infinitely localized solid. The use of two order param- 
eters, average density p and crystallinity x, will allow us 
to explore different pathways from the gas/liquid to the 
solid. Using two order parameters thus provides a richer 
space of possible behaviors and intermediate states than 
does a single order parameter 
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In order to use DFT in practical calculations, it is of 
course necessary to know the functional F [p] . Good ap- 
proximations exist for this functional for the special case 
of hard-sphere interactions but the extension of these to 
other potentials has proven difficult. For this reason, 



where f} is the inverse temperature, F^^ [p; d] is the free 
energy functional for a hard-sphere system with hard- 
sphere diameter d [p] , AF (p) is the difference in free 
energy of a liquid at density p and that of a hard- 
sphere liquid at the same density. In the integral, 
Ac2 (l^i,~?^2; [px] ,d[px]) is the difference in 2-body di- 
rect correlation functions (DCFs) for the interacting sys- 
tem and the hard-sphere system for a density px (T*) = 
p + A (p (T^) — p). The DCFs are not known exactly and 
so it is necessary to introduce approximations to proceed. 
Motivated by the fact that in applications of thermody- 
namic perturbation theory to simple fluids, the correction 
to the hard-sphere free energy is typically similar for FCC 
solids and liquids, we will make the simplest approxima- 
tion which is to assume that the contribution of the third 
term is insignificant. This model has been shown to work 
well for the LJ potential 10] while using the more detailed 
model of Curtin and Ashcroft /gj , which is closely tied to 
the LJ potential, we indeed find the third term to con- 
tribute little. More sophisticated approximations will be 
discussed in a future publication. Here, our interest is not 
the further development of DFT but in its application to 
the question of the kinetics of crystallization. 

In the following, we use the first order WCA per- 
turbation theory[2 0, 0, Q as modified by Ree et 
alE E [13 to calculate the free energy of the liq- 
uid phase. This theory is known to be very accurate 
for a wide class of potentials. The liquid-phase hard- 
sphere diameter calculated from this theory is used for 
both the liquid and the solid phases so that it is in- 
deed solely a function of the average density. For the 
hard-sphere free energy functional we use the fundamen- 
tal measure theory (FMT), specifically the "White Bear 
" functional jisll which gives a good description of the 
hard-sphere phase diagram, in particular reproducing the 
Carnahan-Starling equation of state for the hard-sphere 
liquid. The use of the FMT free energy model is criti- 
cal: previous attempts to perform similar studies made 
use of effective liquid approximations which do not work 
well when the occupancy is treated as a free variable 
and which therefore had to be modified in an ad hoc 



liquid- and solid-state perturbation theory are often used 
as a means of using the hard-sphere theory to approxi- 
mate the functional for other systems 0- Indeed, simple 
manipulations yield the exact relation 



(3) 



manner Ulllillll. The FMT have built into them the 
critical feature that they are sensitive to the local density 
and correctly cause the free energy to diverge if the local 
occupancy grows above one. 

The calculations presented here were performed using 
the standard LJ potential vlj {r) — 4e ^(f)^^ — (f)^^ , 
widely used as a model for the interactions of atomic flu- 
ids, and the potential of tWF VtwF (f) which is intended 
to model the interactions of globular proteins llj . The lat- 
ter consists of a hard-core of diameter a and a modified 
LJ tail VtwF {r) = i'Lj(A^/^(r^ -ct^)^/^) for r > a. where 
A controls the range of the interaction; following ref . , 
we take A = 50 . Figure Q shows the phase diagrams for 
both interaction models calculated using our simplified 
DFT. (Note that the DFT also predicts a spinodal hue 
but, for clarity, it has not been shown.) In both cases, the 
phase diagrams are reproduced surprisingly well given 
the simple models used. The observed deviations from 
the simulation data can be at least partly explained as 
arising from the use of the liquid-state free energy differ- 
ence for the solid which requires knowledge of the liquid 
state at high densities for which even the input hard- 
sphere equation of state is not reliable. Furthermore, the 
determination of phase coexistence is a very sensitive test 
since it depends on getting both the absolute magnitude 
and the slope of the free energies correct. To put this in 
perspective, deviations are observed in Fig^between the 
calculated and simulated gas-liquid coexistence curve for 
the LJ system even though the liquid-state perturbation 
theory gives free energies which differ from simulation by 
less than 1% 

The difference between the phase diagrams of the two 
interaction models is significant and generic. Whereas 
the LJ interaction gives rise to a typical phase diagram 
with a critical point and, at lower temperatures, distinct 
gas and liquid phases and a triple point, the tWF inter- 
action model gives only a single liquidus phase which is 
typical of short ranged interactions. Indeed, as the pa- 
rameter A in the tWF potential is varied from A = 1 to 
A = 50 , the phase diagram evolves continuously from one 



(3F[p] = f3F"'[p-d]+pAF{-p) 



(p(^i) - P) (p(^2) - p) / (1 - A) Ac2 ( 7^1,7^2; [pa] ,rf[pA])rfAdT*2dT*i. 
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FIG. 1: The phase diagram for a LJ potential, left, and the ten 
Wolde-Frenkel potential, right. The full and dashed lines are 
from the model and the points are from simulation, ref. |22ll2^ 
and Q respectively. The dotted lines connect the coexistence 
points used in Figs. 2 and 3. 

similar to the LJ phase diagram to that shown here pos- 
sessing a metastable gas- liquid transition j2^]. This model 
is motivated in part by the fact that a dense metastable 
hquid phase is in fact experimentally observed for some 
proteins. It is this metastable phase which tWF showed 
to play a role in nucleation of the soHd phase from the 
gas. 

Given a reasonable model for the DFT free energy 
functional, we now turn to the question of the effect of 
different paths through density space from a gas of den- 
sity 'Pgi^s and crystallinity Xgas = to a solid with density 
V solid ^'^d crystallinity XsoUd- Here, we consider two can- 
didate paths. The first corresponds to a simultaneous 
densification and ordering of the gas into a solid and is 
parameterized as 

P {x) = Vgas + X {PsoHd - Pgas) (4) 
X{x) = XXsolid 

where x G [0, 1] is an abstract reaction coordinate. This 
might be thought of as the " classical " path. The second 
path we consider is a two step process consisting of first a 
densification at zero crystallinity followed by an ordering 
at fixed density 

P {x) = [Pgas + {PsoHd ~ Pgas)] ® ~ ^) 

+PsoHd'd - 

Xix) = Q (^x - ^ {2x -l)xsoUd- 

Figure |21 shows the free energy landscapes encountered 
using the tWF potential along both paths for coexisting 



FIG. 2: The Gibbs free energy barriers, per atom, for the tWF 
potential. The solid curve is for the classical path and the 
broken curve results from densification followed by ordering. 



gas and solid densities at three different temperatures. 
In all three cases, the classical path requires overcoming 
a free energy barrier as expected. The behavior along 
the non-classical path is more complex. At the highest 
temperature, which lies about the critical point of the 
metastable gas-liquid transition, the non-classical barrier 
is somewhat reduced but the effect is not significant. For 
the intermediate temperature, which is somewhat below 
the critical point, a second, metastable state appears and 
the single free energy barrier splits into two lower barri- 
ers. At the lowest temperature, the barriers encountered 
along the non-classical path are even lower so the ad- 
vantage of this path is even greater. This picture agrees 
well with that developed by Vekilov and co-workers who 
have observed, by means of dynamic light-scattering, the 
presence of short-lived dense liquid droplets in protein 
solutions 'g'I. The results for the LJ system, shown in 
Fig. 121 arc unexpected in that a similar phenomenon 
is observed although the details differ. Again, we show 
three temperatures which are this time all below both the 
critical point and the triple point. The highest temper- 
ature is only just below the triple point and again, it is 
clear that the non-classical path is energetically favored 
relative to the classical path and that this correlates with 
the presence of a metastable dense-liquid state. Unlike 
the previous case, the advantage of passing along the 
non-classical path remains more or less constant as the 
temperature is decreased. Also different is the fact that 
the barrier between the metastable state and the solid 
state is much lower than that between the gas and the 
metastable liquid. This suggests that the droplets in the 
metastable state will crystallize quickly and will be corre- 
spondingly shorter-lived. To check this surprising result, 
we repeated our calculations using the model described in 
ref.[9j,[l9j which is tuned to the LJ potential. While the 
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suits for the free energy landscape for transitions near 
the coexistence lines. For denser gases, which are super- 
saturated with respect to crystallization, the free energy 
of the gas phase moves up so the first barrier is smaller. 
However, it is important to notice that the crossing of the 
free energy barriers is a fundamentally non-equilibrium 
process so that kinematics also plays an important role in 
determining the overall nucleation rates 0,0. Neverthe- 
less, if crystallization kinetics occurs reasonably close to 
equilibrium the free energy functional will play a central 
role since the rates of change of the order parameters will 
be given by the product of its gradient and of a matrix 
of phenomonological parameters. 



FIG. 3: Same as Fig. 2 for the LJ potential. 

barriers were somewhat smaller, the qualitative results 
were the same. 

We have presented calculations of the free energy land- 
scape as a function of density and crystallinity based on 
a simple, robust free energy density functional. Our cal- 
culations involve no input or parameterization except 
for the interaction potential. In both cases studied, 
a model protein and a simple liquid, our results pro- 
vide direct support for the Ostwald rule of stages for 
nucleationji^ since the free energy barriers associated 
with the metastable intermediate states are lower than 
those for a direct transition from gas to solid. The results 
for the model protein interaction agree with the generally 
accepted picture that crystallization proceeds via a two- 
step process of densification followed by crystallization^, 
even when the temperature is slightly above the criti- 
cal point and no metastable intermediate phase exists. 
Interestingly, we find similar behavior for simple fluids 
below the triple point suggesting that crystallization in- 
volving passage through a metastable disordered state 
may be a generic phenomenon. However, in that case, 
the metastable state is expected to be shorter-lived com- 
pared to the nucleation time thus making its experimen- 
tal detection more challenging. The only evidence we 
are aware of for the two-step nucleation mechanism for 
non-protein fluids comes from nonphotochemical laser- 
induced nucleation of small organic molecules [2^ EtI 
and recent molecular dynamics simulations of the crystal- 
lization of AgBr from solution[2^. The short lifetime of 
the metastable phase predicted here would explain why 
it has not so far been observed experimentally in simple 
fluids. 

For the protein model, our results indicate similar bar- 
riers for the gas-liquid and liquid-solid transitions and it 
should be noted that the only experimental results in- 
dicate that the latter should be much higher than the 
formerly. In part, this is because we have presented re- 



A point which could cause concern is that the paths 
through parameter space might cross the spinodal and 
so pass through the two-phase region where it could be 
thought that the use of DFT is problematic. In fact, aside 
from the minima, all points on the curves shown in Fig. 
2 are thermodynamically unstable. However, the funda- 
mental idea underlying DFT is that any density profile 
can be stabilized by means of an external field and that 
the free energies calculated are the intrinsic contribution 
to the free energy when such a stabilizing field is present 
(see, e.g. ref. jljj). Only the intrinsic contribution is 
used here to estimate the barriers as the real system must 
pass through these states without the presence of such a 
stabilizing field. 

One question not answered yet is whether the partic- 
ular pathways discussed here are the optimal, i.e. mini- 
mum energy, pathways. Simple contour plots of the free 
energy show that the non-classical paths used here are in- 
deed very close to the optimal paths as will be discussed 
at length in a future publication. Another question is the 
role of surface tension which should in general increase 
the free energy barriers, as well as the free energies of 
clusters in the metastable and solid states. Since the 
penalty due to surface tension is expected to increase as 
the system passes from the gas to the metastable state to 
the solid state, we expect that the barriers and the free 
energy minima will be shifted accordingly but it seems 
unlikely that the overall picture would change since this 
would require that the addition of surface tension af- 
fect the classical path less than the non-classical path. 
A definitive answer to this question will require calcula- 
tions of free energies for inhomogeneous states which we 
are currently pursuing. 
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